# File name: 	biv_chao_1998_I.R
# By: 		tw
# Purpose:		bayesian IV according to Chao, Philips (1998 JoE) with AdMit package
			version with Om = I
# Last Version:	17.02.09

setwd("D:/Nauka/_R")

library(foreign)
library(fBasics)
library(AdMit)
library(rgl)

d <- read.dta("mondaytruth.dta")

load("biv_chao_1998_I.RData")
load("biv_chao_1998_I_weak.RData")

# ivregress 2sls Y_IV (D_heli = Z_heli), first
# ivregress 2sls Y_IV_weak (D_heli_weak = Z_heli), first
# ivreg Y_IV (D_heli = Z_heli), first
# ivreg Y_IV_weak (D_heli_weak = Z_heli), first
# reg Y_IV_full Z_heli

set.seed(123456)

'KERNEL' <- function(par,Y,D,Z,log=TRUE){
if (is.vector(par)) par <- t(par)

    'KERNEL.SUB' <- function(par,Y,D,Z){

	# logLikelihood function first
		v <- matrix(0,length(Y),2)
		v[,1] <- Y - par[3]*par[1] - par[2] - par[4]*par[1]*Z
		v[,2] <- D - par[3] - par[4]*Z

		Om <- matrix(0,2,2)
		Om[2,1] <- 1
		Om[1,2] <- 1				

		ll <- -0.5*tr(t(v)%*%v)
		
	# log of Prior now

		lp <- .5*log(abs(par[4]^2))
	
		return(ll+lp)
    }
    r <- apply(par,1,KERNEL.SUB,Y=Y,D=D,Z=Z)
    if (!log) r <- exp(r)
    return(r)
}

'G' <- function(par){
	return(par)
}
'H' <- function(par){
	p <- par^2
	return(p)
}

#p <- matrix(0,1000,4)
#p[,1] <- runif(1000,0,.1)
#p[,2] <- runif(1000,15,15.5)
#p[,3] <- runif(1000,0,.2)
#p[,4] <- runif(1000,.5,1)


# First setup d$Y_IV,d$D_heli,d$Z_heli - strong instrument!!!

#k <- KERNEL(p,d$Y_IV,d$D_heli,d$Z_heli)
#print(k[which.max(k)])
#par <- p[which.max(k),]
par <- c(0.02236423, 15.09986261, 0.09221081, 0.84413780)
#par <- c(.047306,15.17836,.102459,.7998847)
#par <- c(0.08166995, 15.08989957,  0.10522812,  0.69945933)

KERNEL(par,d$Y_IV,d$D_heli,d$Z_heli)
#
t1 <- proc.time()
iv.admit <- AdMit(KERNEL,mu0=par,control=list(trace=TRUE,trace.mu=TRUE),Y=d$Y_IV,D=d$D_heli,Z=d$Z_heli)
t2 <- proc.time()

par.rmit <- rMit(N=10000, mit=iv.admit$mit)
#plot(density(par.rmit[,1]))
#x11()
#plot(density(par.rmit[,4]))

t3 <- proc.time()
iv.admitis <- AdMitIS(N=50000,KERNEL,G=G,mit=iv.admit$mit,Y=d$Y_IV,D=d$D_heli,Z=d$Z_heli)
t4 <- proc.time()
print(iv.admitis)

iv2.admitis <- AdMitIS(N=50000,KERNEL,G=H,mit=iv.admit$mit,Y=d$Y_IV,D=d$D_heli,Z=d$Z_heli)
print(iv2.admitis)

t5 <- proc.time()
iv.admitmh <- AdMitMH(N=50000,KERNEL,mit=iv.admit$mit,Y=d$Y_IV,D=d$D_heli,Z=d$Z_heli)
t6 <- proc.time()
print(iv.admitmh)



plot(density(iv.admitmh$draws[,1]))
x11()
plot(density(iv.admitmh$draws[,2]))

t.admit <- t2-t1
t.admitis <- t4-t3
t.admitmh <- t6-t5


#save(iv.admit,par.rmit,iv.admitis,iv2.admitis,iv.admitmh,t.admit,t.admitis,t.admitmh,file="biv_chao_1998_I2.RData")
save(iv.admit,par.rmit,iv.admitmh,t.admit,t.admitmh,file="biv_chao_1998_I_final.RData")



# Second setup d$Y_IV_weak,d$D_heli_weak,d$Z_heli - weak instrument!!!


#k <- KERNEL(p,d$Y_IV_weak,d$D_heli_weak,d$Z_heli)
#print(k[which.max(k)])
#par <- p[which.max(k),]
#par <- c(0.08514918,15.44556057,0.11612739,0.50282896)

#2sls results:
par <- c(.4882475,15.31847,.3012295,.1518955)
#par <- c(.4,15.3,.3,.1)

KERNEL(par,d$Y_IV_weak,d$D_heli_weak,d$Z_heli)

t1 <- proc.time()
wiv.admit <- AdMit(KERNEL,mu0=par,control=list(trace=TRUE,trace.mu=TRUE),Y=d$Y_IV_weak,D=d$D_heli_weak,Z=d$Z_heli)
t2 <- proc.time()

# play(sin(1:10000/20))

parw.rmit <- rMit(N=10000, mit=wiv.admit$mit)
#plot(density(par.rmit[,1]))
#x11()
#plot(density(par.rmit[,4]))

t3 <- proc.time()
wiv.admitis <- AdMitIS(N=50000,KERNEL,G=G,mit=wiv.admit$mit,Y=d$Y_IV_weak,D=d$D_heli_weak,Z=d$Z_heli)
t4 <- proc.time()
print(wiv.admitis)


wiv2.admitis <- AdMitIS(N=50000,KERNEL,G=H,mit=wiv.admit$mit,Y=d$Y_IV_weak,D=d$D_heli_weak,Z=d$Z_heli)
print(wiv2.admitis)

t5 <- proc.time()
wiv.admitmh <- AdMitMH(N=50000,KERNEL,mit=wiv.admit$mit,Y=d$Y_IV_weak,D=d$D_heli_weak,Z=d$Z_heli)
t6 <- proc.time()
print(wiv.admitmh)

plot(density(wiv.admitmh$draws[,1]))
x11()
plot(density(wiv.admitmh$draws[,2]))
#plot.ts(wiv.admitmh$draws[,1])

tw.admit <- t2-t1
tw.admitis <- t4-t3
tw.admitmh <- t6-t5

#save(wiv.admit,parw.rmit,wiv.admitis,wiv2.admitis,wiv.admitmh,tw.admit,tw.admitis,tw.admitmh,file="biv_chao_1998_I_weak1.RData")
save(wiv.admit,parw.rmit,wiv.admitmh,tw.admit,tw.admitmh,file="biv_chao_1998_I_weak_final.RData")








# Robustness check

#d$Y_IV_weak,d$D_heli_weak,d$Z_heli
da <- matrix(0,10000,4)
da[,1] <- d$Y_IV_weak
da[,2] <- d$D_heli_weak
da[,3] <- d$Z_heli
da[,4] <- runif(10000)

q <- da[,4]
dat <- matrix(0,length(q[q<=.1]),3)
j <- 1
for (i in 1:10000){
	if (da[i,4]<=.1) {
		dat[j,] <- da[i,1:3]
		j <- j+1
	}
}

par <- c(.4882475,15.31847,.3012295,.1518955)
KERNEL(par,Y=dat[,1],D=dat[,2],Z=dat[,3])

t1 <- proc.time()
rwiv.admit <- AdMit(KERNEL,mu0=par,control=list(trace=TRUE,trace.mu=TRUE),Y=dat[,1],D=dat[,2],Z=dat[,3])
t2 <- proc.time()

rparw.rmit <- rMit(N=10000, mit=wiv.admit$mit)

t5 <- proc.time()
rwiv.admitmh <- AdMitMH(N=50000,KERNEL,mit=rwiv.admit$mit,Y=d$Y_IV_weak,D=d$D_heli_weak,Z=d$Z_heli)
t6 <- proc.time()
print(wiv.admitmh)
rw.admitmh <- t6-t5

save(dat,rwiv.admit,rparw.rmit,rwiv.admitmh,file="biv_chao_1998_I_weak_robust_final.RData")



